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Abstract 

This article investigates penetrative turbulence in the at¬ 
mospheric boundary layer. Using a large eddy simu¬ 
lation approach, we study characteristics of the mixed 
layer with respect to surface heat flux variations in the 
range from 231.48 Wm“^ to 925.92 Wm~^, and ob¬ 
serve that the surface heterogeneity on a spatial scale 
of 20 km leads to downscale turbulent kinetic energy 
cascade. Coherent fluctuations of mesoscale horizon¬ 
tal wind is observed at 100 m above the ground. Such 
a surface induced temporal oscillations in the horizon¬ 
tal wind suggest a rapid jump in mesocale wind fore¬ 
casts, which is difficult to parameterize using traditional 
one-dimensional ensemble-mean models. Although the 
present work is idealized at a typical scale (20 km) of 
surface heterogeneity, the results help develop effec¬ 
tive subgrid scale parameterization schemes for classical 
weather forecasting mesoscale models. 


1 Introduction 


Earth’s surface is assumed to pla y a distinctiv e influence 
in weather and climate systems (Stub, 19761) . The sur¬ 
face influence generates horizontal convection as a re¬ 
sult of differential heating along one horizontal bound¬ 
ary of the atmospheric boun dary layer (ABL) (ILane . 
2008 : Scotti and Whitel 201 lb . Coherent mesoscale mo¬ 
tions are reported by many studies, where the effect of 
surface heterogeneity are investigated numerically on 
scales of few hundred kilometers, using a grid spac¬ 
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ing of Ax ^ 0(10 km) dSkamarock and Klempl 20081 


Linl [20071). However, in the geophysical fluid dynam¬ 


ics community, controversial opinions exist on whether 
horizontal convection due to differential heating at the 
same level would lead to t urbulence, or drive large scale 
overturning circ u lation s (jPaparella and Youngl 2002 ; 


Scotti and White! 2011). A c entury old experimen¬ 


tal result by ISandstromI (119081) leads to the opinion 
that a sustained circulation cannot occur if the level of 
the heating source is the same as that of the cooling 
source. This hypothesis is supported theoretically by th e 


anti-turb ulence theorem of Paparella and Youngl (l2002h . 


However, Scotti and White ( 201 ih employed a direct nu¬ 
merical simulation of horizontal convection at Rayleigh 
number 10^° to show that a flow driven by the horizontal 
convection exhibits the characteristics of a true turbulent 
flow. Clearly, some fundamental aspects of solar heating 
at the earth’s surface is not fully understood. 


Much of our current understanding of surface induced 
phenomena depends on numerical simulations and mea¬ 
surements/observations. For example, lightning data 
over Houston, USA between years of 1989 and 2000 
indicates that highest flash densities occur over urban 
areas. Measurements of turbulence for the 43 years pe¬ 
riod from 1958 to 2001 in the North Atlantic, USA, and 
European sectors conclude that clear-air turbulence in¬ 
creases in these regions by 40 — 90% as a result of the 
urban induced impact on vertical transport and mixing 
of penetrative ABL turbulence. Note that the average 
horizontal flux of kinetic energy at a height ^ 100 m 
from the ground is about 1000 W m~^. A destabiliza¬ 
tion of this huge energy flux due to differential heating 
on the surface may cause catastrophic impact on the tur- 
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bulent atmosphe re. Since pe netrative turbulence in the 
atmosphere (see LStullLll976h is primarily characterized 
by solar heating, perturbations to boundary layer struc¬ 
tures by human activities (e.g. urba nization) are sensi¬ 
tive to vertical mix ing and transport (IBrvan and Fritsch , 
2002 : Land 2008 ). 

This article reports on a large eddy simulation (LES) 
based numerical model in order to characterize penetra¬ 
tive turbulence over heterogeneous land surface. Pri¬ 
mary objectives of this short article includes the in¬ 
vestigation of evidences whether horizontal convection 
transports large quantities of heat, as well as sustains 
large amounts of diapycnal mixing with a relatively 
small amount of dissipation. We demonstrate that dif¬ 
ferential surface heating on scales of a typical modern 
city cascades kinetic energy downscale, which - in the 
absence of shear - is sufficient to initiate a turbulent 
flow. The computational fluid dynamics approac h used 


in this work is fully detailed bv lAlam et al 


12014. This 


method filters the governing equations using a multires¬ 
olution approach, and employs a Smagorinsky type eddy 
viscosity model for the subfilter scale processes. We use 
a sixth order weighted residual collocation method that 
has no inherent numerical dissipation, and thus, there is 
no need to apply artificial dissipation to ensure numeri¬ 
cal stability. 


2 Methodology 

2.1 Theory 

Penetrative turbulence occurs when an unstably strati¬ 
fied large body of fluid underlies a stably stratified fluid 
layer. In this situation, turbulent eddies are driven ver¬ 
tically by the buoyancy force , and attemp t to penetrate 
into the overlying stable fluid (IStullLll976l) . In the stable 
region, when an eddy reaches its level of buoyancy, it 
returns back downward into the unstable region. How¬ 
ever, due to the gained momentum, eddies often over¬ 
shoot their level of buoyancy, and thus, internal waves 
are excited from the interface between the overlying sta¬ 
bly st ratified fluid and the un derling unstable mixing 
layer ( Stull . 1976 : Lanel 2008h . These waves transports 
kinetic energy to the upper atmo sphere, and may i nitiate 
upper level clear air turbulence. IPeardorff et al.l (Il969h 


investigated such penetrative turbulence experimentally. 
Here, we develop a large eddy simulation model to sim¬ 
ulate the above mentioned phenomena. 

2.2 Equations 

For a compressible atmospheric model, the temperature 
(T) and the pressure (p) are represented by the poten- 
tial temperature (0) and th e Exner function (tt') (e.g.. 


Smolarkiewicz et al.L 12014 . 


0 = T 


(j) 


Rd / Cp 


r = I —pO 
.Po 


Rd/icp — Rd) 


respectively, where po is a reference pressure, Rd is the 
gas constant, Cp is the specific heat at constant pres¬ 
sure, and p is the density. The continuity equation is 
replaced with Eq ([T]i. The notation (xi, cca) = (x, z) and 
(ui, U 3 ) = (m, w) are adopted for simplicity. A spatial 
filter is applied to the momentum and energy equations, 
where Ui{= (ui)) and 6{= (6)) represent filtered ve¬ 
locity and temperature, respectively. More specifically, 
Ui = Uj+u', and0 = 0n+9(z)+O+O' JPeardorff . 197d 
198r)h . Note the separation of the reference temperature 
00 from the background temperature 0{z), which is con¬ 
venient for satisfying the surface condition. The follow¬ 
ing equations are solved in the present work; 


Dtt' fdui 
Dt ^ dxi ’ 


Du^ _ a X 

— CpC/Q „ 'a o J 

Ut OXi UQ OXj 


DO de 

15i ~ ~'^d^ 


dr. 


ej 


dxi 


( 1 ) 


( 2 ) 


(3) 


The simulation region is a vertical plane (x, z) that ex¬ 
tends 100 km horizontally (—50 < x < 50) and 2 km 
vertically (0 < z < 2). A city of scale 20 km exists for 
—10 < X < 10, which is surrounded by rural areas. Re¬ 
cent literature indicates that the most appropriate model 
for the near surface penetrative turbulence is not fully 
unders tood. As a compromise, as discussed by IPope 
(l2000h . we have adopted a resolution that is finer than 
that is used by LES models of the ABL. The finest res¬ 
olution uses Ax = 97.65 m and Az = 3.9 m. Thus, a 
significant fraction of the energy containing large eddies 
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is resolved. Here, some advantages of the costly three- 
dimensional simulation are sacrifi ced for a hig h resolu¬ 
tion two-dimensional idealization (Lane, 20081) . 


turbulence for 6 values of Hq. A purpose of these simu¬ 
lations is to understand the sensitivity of high-amplitude 
surface heat-flux heterogeneity for temporal oscillation 
in mesoscale atmospheric circulations. 


2.3 Large eddy simulation 


The subgrid scale turbulent stress is estimated by 
the popular Smagorinsky (1963) model, = 

—2((7sA)^|5'|S'y, where \S\ = yj2SijSij and Sij = 


(1/2)(|^ + define the filter width by A = 

\/Ax ■ Az, and take Cs = 0.18 for the Smagorinsky 

constant. _ 

According to Deardorfl ( 1980l) . the subgrid scale eddy 
coefficients may be computed by Kh = {1 + 21/A)Km 
and Km = 2((7sA)^|S'|, where Z < A is a subgrid scale 
mixing length. The length scale, I, is related to the sub¬ 
grid scale turbulence energy e' and the buoyancy fre¬ 
quency, N'^ = i-e. I = 0.76\/^/N. If the scale 

of resolved eddies is A = 20 m, then a turbulent Prand- 
tle number Vr = 0.71 gives Z 4 m. The SG S flux for 
buoyancy is rgi = —Kh-§§- ( Deardorfl . 1980h . 


2.4 Surface heat-flux variation 

We consider a time independent profile for the surface 
heat flux variation H^^^,{x) = (Tf^fc) that is given 
by 

^sfc(^) “ (iTsfc)-|-^[tanh^(x-|-A/2)—tanh^(a;—A/2)]. 

Here, (Tf^fc) “ ^ tn (about 37 W m“^) repre¬ 

sents a domain average heat flux for all x in the range 
from —50 to 50 km, where there is no city induced 
heat flux Hq. A is the characteristic wavelength for the 
surface heat flux variation Hq, where the heating re¬ 
gion (city) of the surface is from —10 to 10 km, i.e. 

A = 20 km. This choice for A represents the resolved 
scale for mesoscale numerical weather prediction mod¬ 
els. ^ = 100 is a dimensionless number that leads to 
a continuous sharp interface between the central heat¬ 
ing region and other part of the surface. Thus, the sur¬ 
face heat flux {x) takes approximately the form of 
a square wave without sharp corners, and models the ef¬ 
fect of urban-rural heat flux variation. This article sum¬ 
marizes the sensitivity of surface heat flux on penetrative 


2.5 Validation 


The total heat flux is composed of turbulent (w'9') and 
viscous (ct^) components, where the turbulent compo¬ 
nent vanishes on a flat surface. The dimensionless sur¬ 
face heat flux takes the form — - 


ae\ 


where the 

potential temperature is the same the surface tempera¬ 
ture. 

A scale analysis is necessary to compare the present 
model with that presented by [Pubois and Touzani 
( 2009h . With a fully developed turbulence in the convec¬ 
tive boundary layer, the vertical temper ature profile sat¬ 


isfies ^ = 0 in mixing region (z > 0) (iDeardorff et al. 


I UZ W W « - ' 

19691) . In this case, natural convection heat transfer is 
characterized by the Rayleigh number TZa = ’ 

where H is a vertical length scale, l/^o is the coefficient 
of thermal expansion (1/K), ly is the kinematic viscos¬ 
ity (m^/s), K is thermal conductivity (W/m-K), g is the 
acceleration due to gravity (m/s^), and a is the thermal 
diffusivity (m^/s). In an LES, v and a can be replaced 
with Km and Kh, respectively. So, we set = 10°K 
to match the adiabatic lapse rate (10‘’K/km) of the atmo¬ 
sphere so that ^ = 0. As a result, an increase of by 
a factor of 2 increases TZa by a factor of 10. Thus, the 
values of Hq = 57.87 Wm“^ and 115.74 Wm“^ repre¬ 
sents TZa = 10^ and 10®, respectively. The comparison 
is summarized in Table [T] In addition, we have com¬ 
pared the vertical profile of mean temperature with that 
obtained from the Wangara day 33 experiment, which 
shows an excellent agreement on ^ = 0 between two 
data sets. 


3 Results 

3.1 Surface induced impact of air pollu¬ 
tion 

Fig IHa) demonstrates an example of surface induced 
impact on air pollution. The turbulent plume from the 
lower chimney (75 m tall) moves toward the region of 


3 


































(a) (&) 


Figure 1: (a) Movement of plumes from two chimneys of unequal height as a result of differential surface heating. 
The figure is adapted from Google image database, and represents a pattern of expected circulation, (b) Low level 
converging flow and high level diverging flow is seen from this streamline att = 6.5 h for Jig = 115.74 W m~^. 


warmer surface (city area), and that from the higher 
chimney (150 m tall) moves toward the region of cooler 
surface (urban, ocean). The atmosphere over the warmer 
surface has a decreased thermal stability, which causes a 
low level converging flow. Thus, the air parcels move 
horizontally toward the center of the heated region, 
where they rises upward. The characteristic flow pattern 
in Fig[T]is observed from our numerical simulation. The 
comparison in Fig[T]indicates that our numerical model 
simulates phenomena that is also observed in the nature. 
The streamline plot in Fig[T][6) shows that air moves to 
the left in the lower portion of the boundary layer above 
x > 0, and the air in the upper portion moves to the 
right. This flow pattern is computed at f = 6.5 h from 
the simulation with Hg = 231.48 W m~^. 


3.2 Turbulent mixed layer formation 

To understand the mixing and turbulent transport in 
the absence of shear or mechanically driven turbulence, 
we have initialized the flow with a stable stratification 


Hg 

57.87 Wm-^ 

115.74 Wm"^ 


Present 

D&T 

Present 

D&T 

^min 

-0.0644 

-0.0712 

-0.1672 

-0.1663 

^max 

0.1761 

0.1748 

0.1796 

0.1790 

'^max 

0.2275 

0.2282 

0.3294 

0.3224 

Nu 

0.3262 

0.2951 

0.6892 

0.6435 


Table 1: Com parison of dimensionless e xtreme values 
with that from lDubois and Touzani ( 2009h (D & T). 


and no background wind everywhere except an unstable 
stratification is introduced near the surface. The surface 
Richardson number i?io = controls 

the initial strength of unstable stratification. The flow is 
initialized with Rig = — 1, and allowed to evolve natu¬ 
rally with time. Thus, the resulting cir culations charac 


terize penetrative turbulent convection (iDeardorff et al 


19691) . The time evolution of the streamlines in the en¬ 


tire domain is shown in Fig|2l 

The relative surface heat flux for —10 < x < 10 
is Hg = 925.92 W m’^ (^ 0.75 K m s'^), where the 
background potential temperature 0{z) in the stable re¬ 
gion has a gradient 10”/1 km, and the buoyancy fre¬ 
quency N K, 10“^s“^. As seen from Fig|2l when 
the sun heats the surface, eddies begin to form and rise 
upward. However, they return back downward due to 
stable stratification in the upper atmosphere. This gen¬ 
erates an unstable mixing layer that is adjacent to the 
surface. This turbulent mixed layer underlies a stable 
region aloft. Based on the mean vertical profile of po¬ 
tential temperature (not shown), the depth of the mixed 
layer for this simulation is approximately 800 m. 

When an eddy loses buoyancy, it rises upward, and 
a negative horizontal buoyancy gradient occurs near its 
left edge, which generates a cyclonic circulation; sim¬ 
ilarly, anticycloni c circulation is formed on the other 
edge of the eddy (Lane, 2008t Alaml 201 ih . The pat¬ 
tern of such a circulation and the time evolution of the 
associated span-wise vorticity, during a penetrative tur¬ 
bulent convection, is realized from the streamlines in 
FigEl This vortical pattern play its role as a heat transfer 
agent, thereby causing a imbalance between the buoy- 
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(b) t = 2h 
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Figure 2; The time evolution of span-wise vorticity is demonstrated using the streamline at f = 1 h, f = 2 h, f = 4 h, 
and f = 6 h for the simulation with Hq = 925.92 W m~^. 


ancy force and gravitational force. An important ques¬ 
tion of meteorological interest is whether the process 
leads to horizontal turbulent fluctuation, and if such fluc¬ 
tuation affects local weather prediction. At f = 1 h, 
cyclonic/anticyclonic eddies have been formed near the 
outer edges of the heating region, and have reached a 
height of about 500 m. Later, horizontal convection 
is observed; i.e. turbulent eddies move horizontally, 
where turbulent eddies reach a maximum vertical height 
of about 800 m. Entrainment/detrainment occurs above 
this height. 

3.3 Temporal oscillation and downscale 
energy cascade 

According to the Taylor’s hypothesis, if turbulent statis¬ 
tics is approximately stationary and homogeneous, then 
the turbulent held is advected over the time scales of 
interest. Under this hypothesis, time series of poten¬ 
tial temperature and horizontal velocity, as shown in 
FigS indicate that surface heterogeneity on a scale 
0(20 km) contributes toward downscale energy cas¬ 
cade. The growth of the horizontal potential temperature 
gradient develops a horizontal pressure gradient, which 
in turn generates horizontal wind. The time series of 9 
and u at several vertical locations have been analyzed to 


characterize penetrative turbulence. 

The sensitivity of surface heat flux on the temporal 
fluctuation of 9 is clear from Fig [3 At z = 62.5 m, fre¬ 
quency of oscillation increased at Hq = 925.92 Wm“^ 
compared to Hq — 462.96 Wm“^. However, at 
z = 500 m, turbulence is seen fully developed in both 
cases. This indicates that turbulent kinetic energy cas¬ 
cades downscale as the energy is transported by internal 
waves. To verify that air flow over the simulated city 
is characterized by a horizontal flow of numerous rotat¬ 
ing eddies, we present the horizontal velocity u, which 
is the wind component that is parallel to the direction of 
the surface heat flux variation. The onset of turbulent 
fluctuations is compared between two surfaces fluxes in 

FigEl 


4 Concluding remarks 

In this article, we report some aspects of penetrative 
turbulence in an idealized daytime boundary layer over 
a city that is surrounded by rural areas. As internal 
waves excite from the interface between the mixed layer 
and stable layer, kinetic energy of boundary layer tur¬ 
bulence is transferred to upper label free atmosphere. 
A complete understanding of this mechanism remains 
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Hq = 462.96 Wm-^ ^ = 62.5 m 




Hq = 925.92 Wm-^ z = 62.5 m 



Ho = 925.92 Wm-^ z = 500 m 



Figure 3: Time series of potential temperature 0(0, z, t) near the surface (z = 62.5 m) and near top of the mixing layer 
(z = 500 m). 
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Figure 4: The onset of turbulent fluctuation is shown through time series of horizontal velocity at z = 100 m from the 
surface. 


a challenging future research topic. However, in the 
present article, we have shown that energy cascades 
downscale as turbulence penetrates upward. Our sim¬ 
ulations demonstrate that characteristics of the induced 
horizontal flow may be significantly different depending 
on the surface heterogeneity from the perspective of the 
turbulence that is generated. 

The simulation is idealized. However, we have stud¬ 
ied characteristics of horizontal flow with surface heat- 
flux variation on the scale of 20 km, which is the typ¬ 
ical grid spacing in numerical weather prediction mod¬ 
els. One of our objectives is to understand if thermally 
induced mesoscale perturbation has a connection to pen¬ 
etrative turbulence. In other words, we show that dif¬ 
ferential heating on the same label leads to a sustained 
turbulent flow. Note that we have compromised the 
three-dimensional simulation with a two-dimensional 
one in order to capture the energy containing large ed¬ 
dies. These results suggest that surface heterogeneity 
induced differential heating causes subgrid scale turbu¬ 
lent fluctuations, and the impact of surface heterogene¬ 
ity could be more substantial on mesoscale atmospheric 
flows. In future studies, we plan to continue fully three- 
dimensional LES of penetrative turbulence. We plan to 


investigate how temporal oscillations characterize turbu¬ 
lent energy cascade, particularly with background wind 
conditions and rough surface. 
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